# ---------------------------------------------------------
# TUM - Technichal University of Munich
#
# Original Authors: Taylor Jones
#                   Magdalena Altmann
# Further Authors:  Adrian Wenzel
# Date: 2020-05-06
# Purpose: This source code intends to re-estimate CH4
#          emissions in Munich using a top-down approach.
#          This file provides functions that are used by
#          the main source file.
# Comment: This source code has been developed in 2019,
#          however the date above indicates the first
#          submit to version control.
# ---------------------------------------------------------

require(ggplot2)

# Settings: ------------------------------------------------------------------------------------------------------------------------------
cbPalette <- c("#cc6677","#332288","#ddcc77","#117733","#88ccee","#882255",
               "#44aa99","#999933","#aa4499")

# Quantative Pal:
lumen   <- c('#fefbe9','#f5f3c1','#ddecbf','#c2e3d2','#a8d8dc',
                '#8dcbe4','#7bbce7','#88a5dd','#988ac4','#9a709e',
                '#805770','#46353a')

rainbow <- c('#d1bbd7','#ae76a3','#882e72','#1965B0','#5289C7','#7bafde',
             '#4eb265','#90c987','#cae0ab','#f7f056','#f6c141','#f1932d',
             '#e8601c','#dc050c')


tsj_theme <- theme_bw() + theme(text = element_text(size=15) )
#scale_colour_discrete <- cbPalette
#tsj_theme <- list(tsj_theme, scale_color_manual(values = mycolors) )
theme_set(tsj_theme)

# fig_directory <- 'figs/' ## Depeached defined in Run_Inversion.r

addLegend_decreasing <- function (map, position = c("topright", "bottomright", "bottomleft", 
                                                    "topleft"), pal, values, na.label = "NA", bins = 7, colors, 
                                  opacity = 1, labels = NULL, labFormat = labelFormat(), 
                                  title = NULL, className = "info legend", layerId = NULL, 
                                  group = NULL, data = getMapData(map), decreasing = FALSE) {
  position <- match.arg(position)
  type <- "unknown"
  na.color <- NULL
  extra <- NULL
  if (!missing(pal)) {
    if (!missing(colors)) 
      stop("You must provide either 'pal' or 'colors' (not both)")
    if (missing(title) && inherits(values, "formula")) 
      title <- deparse(values[[2]])
    values <- evalFormula(values, data)
    type <- attr(pal, "colorType", exact = TRUE)
    args <- attr(pal, "colorArgs", exact = TRUE)
    na.color <- args$na.color
    if (!is.null(na.color) && col2rgb(na.color, alpha = TRUE)[[4]] == 
        0) {
      na.color <- NULL
    }
    if (type != "numeric" && !missing(bins)) 
      warning("'bins' is ignored because the palette type is not numeric")
    if (type == "numeric") {
      cuts <- if (length(bins) == 1) 
        pretty(values, bins)
      else bins	
      
      if (length(bins) > 2) 
        if (!all(abs(diff(bins, differences = 2)) <= 
                 sqrt(.Machine$double.eps))) 
          stop("The vector of breaks 'bins' must be equally spaced")
      n <- length(cuts)
      r <- range(values, na.rm = TRUE)
      cuts <- cuts[cuts >= r[1] & cuts <= r[2]]
      n <- length(cuts)
      p <- (cuts - r[1])/(r[2] - r[1])
      extra <- list(p_1 = p[1], p_n = p[n])
      p <- c("", paste0(100 * p, "%"), "")
      if (decreasing == TRUE){
        colors <- pal(rev(c(r[1], cuts, r[2])))
        labels <- rev(labFormat(type = "numeric", cuts))
      }else{
        colors <- pal(c(r[1], cuts, r[2]))
        labels <- rev(labFormat(type = "numeric", cuts))
      }
      colors <- paste(colors, p, sep = " ", collapse = ", ")
      
    }
    else if (type == "bin") {
      cuts <- args$bins
      n <- length(cuts)
      mids <- (cuts[-1] + cuts[-n])/2
      if (decreasing == TRUE){
        colors <- pal(rev(mids))
        labels <- rev(labFormat(type = "bin", cuts))
      }else{
        colors <- pal(mids)
        labels <- labFormat(type = "bin", cuts)
      }
      
    }
    else if (type == "quantile") {
      p <- args$probs
      n <- length(p)
      cuts <- quantile(values, probs = p, na.rm = TRUE)
      mids <- quantile(values, probs = (p[-1] + p[-n])/2, 
                       na.rm = TRUE)
      if (decreasing == TRUE){
        colors <- pal(rev(mids))
        labels <- rev(labFormat(type = "quantile", cuts, p))
      }else{
        colors <- pal(mids)
        labels <- labFormat(type = "quantile", cuts, p)
      }
    }
    else if (type == "factor") {
      v <- sort(unique(na.omit(values)))
      colors <- pal(v)
      labels <- labFormat(type = "factor", v)
      if (decreasing == TRUE){
        colors <- pal(rev(v))
        labels <- rev(labFormat(type = "factor", v))
      }else{
        colors <- pal(v)
        labels <- labFormat(type = "factor", v)
      }
    }
    else stop("Palette function not supported")
    if (!any(is.na(values))) 
      na.color <- NULL
  }
  else {
    if (length(colors) != length(labels)) 
      stop("'colors' and 'labels' must be of the same length")
  }
  legend <- list(colors = I(unname(colors)), labels = I(unname(labels)), 
                 na_color = na.color, na_label = na.label, opacity = opacity, 
                 position = position, type = type, title = title, extra = extra, 
                 layerId = layerId, className = className, group = group)
  invokeMethod(map, data, "addLegend", legend)
}


# Plot the em27 observations ---------------------------------------------------------------------------------------------------------
plt_obs <- function( y_df , save_plot = T , fig_directory. = fig_directory){
  output <- ggplot(y_df,aes(x=recep,y=ob,color=em)) +
    theme_gray() +
    geom_point() + 
    geom_line() + 
    facet_wrap(~date,scales="free_x") +
    labs(x="Observation Time [UTC]",y="XCH4  [ppm]",color="Instrument") +
    theme(text = element_text(size=15) ) +
    scale_colour_manual(values=cbPalette, breaks=ems)
  if(save_plot){
    #ggsave( paste0( fig_directory , '/observations_', y_df$date[1], '.png') )
    ggsave( paste0( fig_directory,'/observations_', date, '.png') )
  }
  return(print(output))
}


# Plot the footprints for each em and sector, all stacked up --------------------------------------------------------------------------
plt_stacked <- function( y_df , save_plot = T , fig_directory. = fig_directory ){
  output <- ggplot(y_df,aes(x=recep,y=fp*1000,fill=sector)) + 
    theme_gray() +
    theme(text = element_text(size=15) ) +
    geom_area() + 
    facet_wrap(~em,scales="fixed") +
    scale_fill_manual(values=cbPalette, breaks=sectors ) +
    #ylim(0,5)+
    labs(x="Receptor Time [UTC]", y="XCH4 Contribution  [ppb]" ,fill="Sector",title = paste0( "Prior Expecting Contributions: ", date) ) 
    #facet_wrap(~date,scales="free")
  if(save_plot){
    #ggsave( paste0( fig_directory , '/stacked_prior_', date, '.png') )
    ggsave( paste0( fig_directory,'/stacked_prior_', date, '.png') )
  }
  return(print(output))
}

# Plot Emissions bar graph -----------------------------------------------------------------------------------------------------------
plt_bars <- function( bars_df , save_plot = T ){
  output <- ggplot(bars_df, aes(x=type,y=emission)) +
    theme_gray() + 
    geom_bar(position="dodge", stat="identity",width=.5,fill=cbPalette[1]) +
    geom_errorbar(position="dodge", aes(ymin=lower, ymax=upper),width=.2) +
    theme(text = element_text(size=15) ) +
    labs(x="Date",y=bquote("Total CH4 Emissions  ["~Gg~yr^{-1}~"]"),title="") +
    geom_hline(yintercept = prior_totals[1], linetype=2) +
    geom_hline(yintercept = 0) +
    #ylim(0,85)
    #scale_y_continuous(breaks=seq(-40,90,20))
  if(save_plot){
    ggsave( paste0( fig_directory,'/bars_', date, '.png') )
  }
  return(print(output))
}

#Plot the main results: ---------------------------------------------------------------------------------------------------------------
plt_result <- function(df,bk_df, save_plot = T, fig_directory. = fig_directory){
  output <- ggplot(df,aes(x=recep)) +
    theme_gray() + 
    geom_point(aes(y=ob,color=em)) +
    geom_line(aes(y=y_hat,color=em)) +
    #geom_line(aes(y=y_bk,color=em),linetype=2) +
    scale_colour_manual(values=cbPalette, breaks=ems ) +
    geom_line(data=bk_df, aes(x=t,y=b_hat),color="black",linetype=2) +
    scale_x_datetime(labels = date_format("%H"), 
                     breaks = seq( as.POSIXct("2016-05-13"),as.POSIXct("2016-05-14"),by="2 hours" ),
                     limits = c(as.POSIXct("2016-05-13 10:00", tz="UTC"),as.POSIXct("2016-05-13 23:00", tz="UTC") ) ) +
    
    ylim(1.84, 1.86)+
    labs(x="Receptor Time [UTC]", y=bquote(XCH[4]~"[ppm]") ,color="Sensor ID") +
    #ggtitle("observations (points) + opt. total column values (solid line) \n + opt. bkgd at station (dashed line)") +
    theme(text = element_text(size=15) )
    #facet_wrap(~date,scales="free")
  if(save_plot){
    ggsave( paste0( fig_directory,'/result', date, '.png') )
  }
  return(print(output))
}

#Plot the correlation between model and obs: -----------------------------------------------------------------------------------------
plt_corr <- function( df, save_plot = T){
  output <- ggplot(df) + 
    theme_gray() +
    geom_point(aes(x=ob,y=y_hat,color=as.factor(date))) +
    scale_colour_manual(values=cbPalette ) +
    labs(x="Observations [ppm]",y="Model  [ppm]" ,color="Date",title = "") + #bquote(R^2~ '='~ .(round(r2_concentration, 3)) )) +
    #labs(x="Observations [ppm]",y="Model  [ppm]" ,color="Date") +
    theme(text = element_text(size=15) ) +
    coord_equal() +
    geom_abline( slope = 1,intercept = 0)
  if(save_plot){
    ggsave( paste0( fig_directory,'/corr_', date, '.png') )
  }
  return(print(output))
}

# Plot the correlation between model and obs, but first subtract the model background from both: -------------------------------------
plt_corr_no_bkgd <- function( df , save_plot = T){ 
  output <- ggplot(df) + 
    theme_gray() +
    geom_point(aes(x=1000*(ob-y_bk),y=1000*(y_hat-y_bk),color=as.factor(date))) +
    scale_colour_manual(values=cbPalette ) +
    labs(x="Observations [ppb]",y="Model  [ppb]" ,color= "Date",title = "" ) +# bquote(R^2~ '='~ .(round(r2_enhancement, 3)) )) +
    #labs(x="Observations [ppb]",y="Model  [ppb]" ,color="Date") +
    theme(text = element_text(size=15) ) +
    coord_equal() +
    geom_abline( slope = 1,intercept = 0) +
    xlim(-4,4) +
    ylim(0,3)
  if(save_plot){
    ggsave( paste0( fig_directory,'/corr_no_bkgd_', date, '.png') )
  }
  return(print(output))
}

# Plot just the backgrounds: ----------------------------------------------------------------------------------------------------------
plt_bkgd <- function(df,bk_df, save_plot = T, fig_directory. = fig_directory){
  output <- ggplot(df,aes(x=recep)) +
    theme_gray() +
    #geom_point(aes(y=y_bk,color=em)) +
    geom_line(aes(y=y_bk,color=em)) +
    scale_colour_manual(values=cbPalette, breaks=ems ) +
    #geom_point(data=bk_df, aes(x=t,y=b_hat),color="black") +
    geom_line(data=bk_df, aes(x=t,y=b_hat),color="black",linetype=3) +
    #ylim(1.83, 1.87)+
    labs(x="Receptor Time [UTC]", y="background XCH4   [ppm]" ,color="Instrument") +
    #ggtitle("optimized bkgd (black, dashed line) \n + optimized bkgd at station (color instr., solid line)") +
    theme(text = element_text(size=15) ) +
    facet_wrap(~date,scales="free")
  if(save_plot){
    ggsave( paste0( fig_directory,'/bkgd_', date, '.png') )
  }
  return(print(output))
}


# Plot modeled + observed + station bkgd: ----------------------------------------------------------------------------------------------------------
plt_bkgd_with_stat <- function(df,bk_df, y_df, b_df){
  output <- ggplot(df,aes(x=recep)) +
    theme_gray() +
    #geom_point(aes(y=y_bk,color=em)) +
    geom_line(aes(y=y_bk,color=em)) +
    scale_colour_manual(values=cbPalette, breaks=ems ) +
    #geom_point(data=bk_df, aes(x=t,y=b_hat),color="black") +
    geom_line(data=bk_df, aes(x=t,y=b_hat),color="black",linetype=2) +
    geom_point(data=y_df, aes(x=recep,y=ob, color = "bkgd station")) +
    geom_line(data=y_df, aes(x=recep,y=ob, color = "bkgd station")) +
    scale_colour_manual(values=cbPalette) +
    #geom_line(data=b_df, aes(x=t,y=b_hat),color="red") +
    labs(x="Receptor Time [UTC]", y="background XCH4   [ppm]" ,color="Instrument") +
    #labs(x="Receptor Time [UTC]", y="background XCH4   [ppm]",color="") +
    theme(text = element_text(size=15) ) +
    facet_wrap(~date,scales="free")
  return(print(output))
}

# Plot the modeled and observed bkgd at bkgd station: ----------------------------------------------------------------------------------------------------------
plt_bkgd_mod_obs <- function(df,bk_df, y_df){
  output <- ggplot(df,aes(x=recep)) +
    theme_gray() +
    #geom_point(data=bk_df, aes(x=t,y=b_hat,color="opt. bkgd")) +
    geom_line(data=bk_df, aes(x=t,y=b_hat),color = "black",linetype=2) +
    geom_point(data=y_df, aes(x=recep,y=ob,color="observations at \nbkgd station")) +
    geom_line(data=y_df, aes(x=recep,y=ob,color="observations at \nbkgd station")) +
    scale_colour_manual(values=cbPalette)+
    labs(x="Receptor Time [UTC]", y="background XCH4   [ppm]",color="") +
    theme(text = element_text(size=15) ) +
    facet_wrap(~date,scales="free")
  return(print(output))
}

# Plot background posterior covariance: -----------------------------------------------------------------------------------------------
plt_bk_pos <- function(bk_df,log=F, save_plot = T){
  require(scales)
  output <- ggplot( bk_df , aes(x=t,y=pos_bk_sigma*1000, color="Posterior") ) +
    theme_gray() +
    theme(text = element_text(size=15) ) +
    geom_line(size=1) +
    facet_wrap(~date, scales = "free_x") +
    scale_x_datetime(labels = date_format("%H:00")) +
    geom_line( aes(x=t,y=prior_bk_sigma*1000, color = "Prior",),size=1 ) +
    labs(x="Background Time [UTC]", y="Model Background Standard Error [ppb]",color="") +
    theme(legend.position = c(0.75, 0.25))
  if(log){
    output <- output + scale_y_continuous(trans = 'log10')
  }
  if(save_plot){
    ggsave( paste0( fig_directory,'/bkgd_pos_covariance_', date, '.png') )
  }
  return(print(output))
}

# Plot Background Error: --------------------------------------------------------------------------------------------------------------
plt_bk_error <- function(y_df, save_plot = T){
  output <- ggplot(y_df, aes(x=recep,y=bk_err*1000,color=em)) +
    theme_gray() + 
    theme(text = element_text(size=15) ) +
    geom_point() +
    geom_line() +
    scale_colour_manual(values=cbPalette ) +
    facet_wrap(~date, scales = "free_x") +
    labs(x="Observation Time [UTC]", y="Model Background Standard Error [ppb]",color="EM") +
  if(save_plot){
    ggsave( paste0( fig_directory,'/bkgd_error_', date, '.png') )
  }
  return(print(output))
}

# Plot total result: -----------------------------------------------------------------------------------------------------------------
plt_total_result <- function(df, save_plot = T){
  output <- ggplot(df,aes(x=recep)) +
    geom_line(aes(y = y_hat*1000, color = 'Inversion'), linetype = 'solid', size = 1) +
    geom_line(aes(y = ob*1000, color = 'Observation'), linetype = 'solid', size = 1) +
    geom_line(aes(y = y_bk*1000, color = 'Background'), linetype = 'dashed') +
    theme_bw() +  theme(text = element_text(size=15)) +
    facet_grid(em ~ date, scales = "free_x") +
    labs(x="Observation Time [UTC]", y="CH4 conc. values [ppb]",color="Legend")
  if(save_plot){
    ggsave( paste0( fig_directory,'/concentration_values_', date, '.png') )
  }
  return(print(output))
}

# Plot enhancement: ------------------------------------------------------------------------------------------------------------------
plt_enhancement <- function(df,bk_df, save_plot = F){
  output <- ggplot(m$df,aes(x=recep)) +
    geom_line(aes(y=(y_hat-y_bk)*1000, color = 'inversion result'), linetype = 'solid', size = 1) +
    #geom_line(aes(y=(y_prior-bkgd_prior)*1000, color = 'prior expect. contribution'), linetype = 'dashed') +
    theme_bw() +  theme(text = element_text(size=15)) +
    facet_grid(em ~ date, scales = "free_x") +
    labs(x="Observation Time [UTC]", y="XCH4 enhancement values [ppb]",color="Legend") 
  if(save_plot){
    ggsave( paste0( fig_directory,'enhancement_values_', date, '.png') )
  }
  return(print(output))
}

plt_y_hat <- function(m){
  output <- ggplot(m) + geom_line(aes(x=recep,y=y_hat*1000,color=em)) +
    geom_line(aes(x = recep,y=y_bk*1000,color = em)) +
    theme_gray() + 
    theme(text = element_text(size=15) ) +
    geom_point() +
    scale_colour_manual(values=cbPalette ) +
    facet_wrap(~date, scales = "free_x") +
    labs(x="Observation Time [UTC]", y="Model values [ppb]",color="EM")
  return(print(output))
}

plt_BIM <- function(bim_df, pltlog=T, save_plot=T , fig_directory.=fig_directory){
  if (pltlog){
    output <- ggplot(bim_df) +
      geom_tile(aes(x=recep,y=bkgd,fill=log(value))) +
      # scale_fill_gradientn(colors=lumen, limits=c(-10,-0.05)) +
      scale_fill_gradientn(colors=lumen) +
      geom_abline( slope = 1,intercept = 0,size=.5) +
      geom_abline( slope = 1,intercept = -60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -2*60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -3*60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -4*60*60,linetype=2,size=.5) +
      scale_x_datetime(labels = date_format("%H"), 
                       breaks = seq( as.POSIXct("2016-05-13"),as.POSIXct("2016-05-14"),by="2 hours" ),
                       limits = c(as.POSIXct("2016-05-13 10:00", tz="UTC"),as.POSIXct("2016-05-13 23:00", tz="UTC") ) ) +
      scale_y_datetime(labels = date_format("%H"), 
                       breaks = seq( as.POSIXct("2016-05-13"),as.POSIXct("2016-05-14"),by="2 hours" ),
                       limits = c(as.POSIXct("2016-05-13 06:00", tz="UTC"),as.POSIXct("2016-05-13 23:00", tz="UTC") ) ) +
      coord_equal() +
      theme_gray() +
      #facet_wrap(~em) +
      labs(x="Receptor Time [UTC]",y="Background Time [UTC]", fill=bquote(log[10]~"(Background Influence)")) +
      theme(text = element_text(size=15) )
  } else {
    output <- ggplot(bim_df) +
      geom_tile(aes(x=recep,y=bkgd,fill=value)) +
      # scale_fill_gradientn(colors=lumen, limits=c(0.00001,0.9)) +
      scale_fill_gradientn(colors=lumen, limits=c(.00001,.4)) +
      geom_abline( slope = 1,intercept = 0,size=.5) +
      geom_abline( slope = 1,intercept = -60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -2*60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -3*60*60,linetype=2,size=.5) +
      geom_abline( slope = 1,intercept = -4*60*60,linetype=2,size=.5) +
      coord_equal() +
      theme_gray() +
      facet_wrap(~em) +
      labs(x="Receptor Time [UTC]",y="Background Time [UTC]",fill="Background Influence") +
      theme(text = element_text(size=15) ) +
      scale_x_datetime(labels = date_format("%H"), 
                       breaks = seq( as.POSIXct("2016-05-13"),as.POSIXct("2016-05-14"),by="4 hours" ),
                       limits = c(as.POSIXct("2016-05-13 10:00", tz="UTC"),as.POSIXct("2016-05-13 23:00", tz="UTC") ) ) +
      scale_y_datetime(labels = date_format("%H"), 
                       breaks = seq( as.POSIXct("2016-05-13"),as.POSIXct("2016-05-14"),by="4 hours" ),
                       limits = c(as.POSIXct("2016-05-13 06:00", tz="UTC"),as.POSIXct("2016-05-13 23:00", tz="UTC") ) )
    
      
  }
  if(save_plot){
    ggsave( paste0( fig_directory,'/BIM_', date, '.png') )
  }
  return(print(output))
}


#different bkgd prior plotting----------------------------------------------------------------

cbPalette <- c("#cc6677","#332288","#ddcc77","#117733","#88ccee","#882255",
               "#44aa99","#999933","#aa4499")

plt_b <- function(bkgd, b_84){
  output <- ggplot(bkgd,aes(x=time)) +
    theme_gray() +
    #geom_point(aes(y=b_hat*1000-b_84*1000,color=bkgd_prior)) +
    geom_point(aes(y=b_hat*1000,color=bkgd_prior)) +
    #geom_line(aes(y=b_hat*1000-b_84*1000,color=bkgd_prior)) +
    geom_line(aes(y=b_hat*1000,color=bkgd_prior)) +
    scale_colour_manual(values=cbPalette) +
    #labs(x="Time [UTC]", y="background difference [ppb]" ,color="bkgd prior value [ppm]") +
    labs(x="Time [UTC]", y="background XCH4 [ppb]" ,color="bkgd prior value [ppm]") +
    theme(text = element_text(size=15) ) 
  return(print(output))
}



# Make a map of all of the priors, as well as the em27 locations and the point sources: ----------------------------------------------
plt_inventory_map <- function(ras,ems,em_lats,em_lons) {
  library(leaflet)
  library(leaflet.opacity)
  library(markdown)
  library(geojsonio)

  fp_color <- c('#eff28c','#95D840FF','#55C667FF','#29AF7FFF','#1F968BFF','#287D8EFF','#33638DFF','#404788FF','#482677FF','#440154FF')

  fp_scale <- c(0.04, 5.7) #adjust to min/max value of the footprint
  #log_scale <- c(-2.5, -0.5) #scale for inventory without pp
  log_scale <- c(-2.5, 0.5) #scale for inventory with pp
  
  pltfp <- colorNumeric(palette=fp_color, domain=fp_scale,na.color = "transparent")
  plt  <- colorNumeric(palette=lumen, domain=log_scale,na.color=lumen[1])


  ras_ex <- extent(ras)
  
  ###load the shape of munich
  munich <- geojsonio::geojson_read(paste0(path_campaign,"/",path_geojson))
  
  
  plt_map <- leaflet() %>% 
    
  #addProviderTiles('Esri.WorldImagery') %>% setView(11.6,48.1, zoom = 9)             #satellite picture
  addProviderTiles(providers$Esri.WorldTopoMap) %>% setView(11.6,48.1, zoom = 9)      #providers$Stamen.TonerBackground, 'Esri.WorldImagery'
  
  #add inventory map
  plt_map <- addRasterImage(plt_map,log10(ras_stack[[1]]), opacity = .75,colors = plt, layerId = "ch4", group = "prior")
  #plt_map <- addOpacitySlider(plt_map, layerId = "ch4")
  plt_map <- leaflet::addLegend(plt_map, pal=plt, values=log_scale, title = "CH4<br>[log(umol/m2s)]")
  
  #add footprints
  #plt_map <- addRasterImage(plt_map,ras_fp*1000, opacity = 0.9,colors = pltfp, layerId = "fp", group = "fp")
  #plt_map <- leaflet::addLegend(plt_map, pal=pltfp, values=fp_scale, title = "footprint <br> [ppb/(umol m-2 s-1]")
 
  #add HKW Sued (power plant) and WWTP location
  #plt_map <- addCircleMarkers(plt_map, 11.55583, 48.11444, radius = 3.5, weight = 1, color = 'black', fillColor = 'blue', label = "PP", labelOptions = labelOptions(noHide = T, textsize = "15px"), stroke=TRUE,fillOpacity = 1, group = "point_source")
  #plt_map <- addCircleMarkers(plt_map, 11.6293331, 48.2088205, radius = 3.5, weight = 1, color = 'black', fillColor = 'blue', label = "WWTP", labelOptions = labelOptions(noHide = T, textsize = "15px"), stroke=TRUE,fillOpacity = 1, group = "point_source")

  #add em locations
  plt_map <- addCircleMarkers(plt_map,em_lons,em_lats,radius = 3.5, weight = 1, color = 'black', fillColor = '#c50032', label = ems, labelOptions = labelOptions(noHide = T, textsize = "15px"), stroke=TRUE,fillOpacity = 1,group = "EM Sites")
  
  #add boundary of Munich
  plt_map <- addGeoJSON(plt_map, munich, group = "Munich boundary", weight = 2, fill = F)
  
  #add boundary of domain 
  plt_map <- addRectangles(plt_map, ras_ex[1], ras_ex[3], ras_ex[2], ras_ex[4], fillColor = "transparent",color = "#444444" )
  #plt_map <- addRectangles(plt_map, 11.17, 47.95, 11.97, 48.35, fillColor = "transparent", color = "#444444") 
  
  plt_map <- addScaleBar(plt_map)
  #plt_map <- addLayersControl(plt_map,overlayGroups = c("prior","EM Sites","HKW", "Munich boundary", "fp"),options = layersControlOptions(collapsed = FALSE))
  
  return(print(plt_map))
}

  

#bs_sf <- bs_sf*prior_total
#bs_df <- data.frame(scaling_factor=bs_sf)
#plt_bs <- ggplot(bs_df,aes(x=scaling_factor)) + 
#  geom_histogram(aes(y=..density..),binwidth=.05,color="black", fill="lightblue") +
#  geom_vline(aes(xintercept=median(scaling_factor)), color="blue", linetype="dashed", size=1) +
#  theme(text = element_text(size=15) ) +
#  labs(x="Scaling Factor",y="Density",title=date)
#p#rint( quantile(bs_sf_1,c(.05,.95)) )
#ba#rs <- rbind( bars , data.frame(date = date, value=median(bs_sf_1),lower = quantile(bs_sf_1,c(.05,.95))[[1]], upper = quantile(bs_sf_1,c(.05,.95))[[2]], sector = sectors_to_bootstrap[1] ))
#

